% Given parameters
C1 = 0.22e-6; % F
C2 = 0.1e-6;  % F
R = 1e3;      % Ohm
L = 10e-6;    % H

% Time span for the solution
tspan = [0, 0.03]; % Adjust the time span as needed

% Initial conditions for V2, V3, V4, V5, integral_V2, integral_V3, integral_V4, integral_V5
V0 = [0, 0, 0, 0, 0, 0, 0, 0];

% Solve the system of ODEs using ode45
[t, V] = ode45(@(t, V) circuitEquations(t, V, C1, C2, R, L), tspan, V0);

% Calculate V1 as a numeric array based on the time array 't'
V1 = sin(2 * pi * 89000 * t); % Now V1 is directly numeric, matching the length of t

% Plot the results in separate figures
figure;
plot(t, V1, 'r');
xlabel('Time (s)');
ylabel('Voltage (V)');
title('V1 vs Time');
grid on;

figure;
plot(t, V(:, 1), 'b');
xlabel('Time (s)');
ylabel('Voltage (V)');
title('V2 vs Time');
grid on;

figure;
plot(t, V(:, 2), 'g');
xlabel('Time (s)');
ylabel('Voltage (V)');
title('V3 vs Time');
grid on;

figure;
plot(t, V(:, 3), 'k');
xlabel('Time (s)');
ylabel('Voltage (V)');
title('V4 vs Time');
grid on;

figure;
plot(t, V(:, 4), 'm');
xlabel('Time (s)');
ylabel('Voltage (V)');
title('V5 vs Time');
grid on;

% Function that defines the system of ODEs
function dVdt = circuitEquations(t, V, C1, C2, R, L)
    V1 = sin(2 * pi * 89000 * t); % V1 is calculated numerically for each time 't'
    V2 = V(1);
    V3 = V(2);
    V4 = V(3);
    V5 = V(4);
    integral_V2 = V(5);
    integral_V3 = V(6);
    integral_V4 = V(7);
    integral_V5 = V(8);

    % Define the differential equations based on the given equations
    dV1_dt = 2 * pi * 89000 * cos(2 * pi * 89000 * t);
    dV2_dt =  -(L*V2 - L*V1 + R*integral_V2 + R*integral_V3 + R*integral_V4 + R*integral_V5 - 178000*C2*L*R*pi*cos(178000*pi*t))/(C2*L*R);
    dV3_dt =  -(C1*L*V2 - C1*L*V1 - C2*L*V2 + C2*L*V3 + C1*R*integral_V2 + C1*R*integral_V3 + C1*R*integral_V4 + C2*R*integral_V3 + C1*R*integral_V5 + C2*R*integral_V4 + C2*R*integral_V5 - 178000*C1*C2*L*R*pi*cos(178000*pi*t))/(C1*C2*L*R);
    dV4_dt = -(C1*L*V2 - C1*L*V1 - C1*L*V3 - C2*L*V2 + C1*L*V4 + C2*L*V3 + C1*R*integral_V2 + C1*R*integral_V3 + 2*C1*R*integral_V4 + C2*R*integral_V3 + 2*C1*R*integral_V5 + C2*R*integral_V4 + C2*R*integral_V5 - 178000*C1*C2*L*R*pi*cos(178000*pi*t))/(C1*C2*L*R);
    dV5_dt = -(C1*L*V2 - C1*L*V1 - C1*L*V3 - C2*L*V2 + C1*L*V4 + C2*L*V3 - C2*L*V4 + C2*L*V5 + C1*R*integral_V2 + C1*R*integral_V3 + 2*C1*R*integral_V4 + C2*R*integral_V3 + 2*C1*R*integral_V5 + C2*R*integral_V4 + 2*C2*R*integral_V5 - 178000*C1*C2*L*R*pi*cos(178000*pi*t))/(C1*C2*L*R);

    % Define the differential equations for the integrals
    dintegral_V2_dt = V2;
    dintegral_V3_dt = V3;
    dintegral_V4_dt = V4;
    dintegral_V5_dt = V5;

    % Return the derivatives as a column vector
    dVdt = [dV2_dt; dV3_dt; dV4_dt; dV5_dt; dintegral_V2_dt; dintegral_V3_dt; dintegral_V4_dt; dintegral_V5_dt];
end



